model {
	for (i in 1:N){
		for (t in 1:7){y[i,t] ~ dinterval(
			z[i,t],cut[1:2])
			e[i,t] <- z[i,t]-h[i,t]
			p2[i,t] <- pow(e[i,t],2)
			znew[i,t] ~ dnorm(h[i,t],1)
			ynew[i,t] <- (step(cut[1]-znew[i,t]) + 	
				2*step(znew[i,t]-cut[1])*step(cut[2]-znew[i,t]) +
				3*step(znew[i,t]-cut[2]))-1
			Match[i,t] <- equals(ynew[i,t],y[i,t])
			for (c in 1:3){
				d[i,t,c] <- equals(y[i,t],(c-1))
				dnew[i,t,c] <- equals(ynew[i,t],(c-1))
			}
		}
			
		for (c in 1:3){
			dsum[i,c] <- sum(d[i,1:7,c])
			dnewsum[i,c] <- sum(dnew[i,1:7,c])
		}

		
		P1[i] <- sum(p1[i,2:7])
		P2[i] <- sum(p2[i,2:7])
					
		z[i,1] ~ dnorm(h[i,1],1)
		
		h[i,1] <- xi[i]*scale +
			delta[1]*age[i,1] + 
			delta[2]*nkids[i,1] + 
			delta[3]*hhsize[i,1] + 
			delta[4]*ownocc[i,1] + 
			delta[5]*hsequi[i,1] + 
			delta[6]*incshr[i,1] + 
			delta[7]*perminc[i,1] + 
			delta[8]*transinc[i,1] + 
			delta[9]*div[i,1] + 
			delta[10]*uempl[i,1] + 
			delta[11]*fem[i] +
			delta[12]*nw[i] +
			delta[13]*degree[i] +
			delta[14]*alvl[i] +
			delta[15]*olvl[i] +
			delta[16]*union[i,1]+
			delta[17]*zg1[i] +
			delta[18]*zg2[i] +
			delta[19]*zg3[i] +
			delta[20]*zg4[i] +
			delta[21]*lond[i]
			
		for (t in 2:7){
			z[i,t] ~ dnorm(h[i,t],1)
			h[i,t] <- beta0 + gamma*z[i,t-1] + xi[i] +
				beta[1]*age[i,t] + 
				beta[2]*nkids[i,t] + 
				beta[3]*hhsize[i,t] + 
				beta[4]*ownocc[i,t] + 
				beta[5]*hsequi[i,t] + 
				beta[6]*incshr[i,t] + 
				beta[7]*perminc[i,t] + 
				beta[8]*transinc[i,t] + 
				beta[9]*div[i,t] + 
				beta[10]*uempl[i,t] + 
				beta[11]*fem[i] +
				beta[12]*nw[i] +
				beta[13]*degree[i] +
				beta[14]*alvl[i] +
				beta[15]*olvl[i] +
				beta[16]*union[i,t]
				
				p1[i,t] <- e[i,t]*e[i,t-1]	
		}
		
		# subgroup indicator
  		xi[i] <- phi[S[i]]
  		S[i] ~ dcat(p[1:M])   
		for (m in 1:M) {memb[i,m] <- equals(S[i],m)}
	}

	for (c in 1:3){
		ppp.obs[c] <- mean(dsum[,c])
		ppp.pred[c] <- mean(dnewsum[,c])
	}

	# base measure G0 
	for (m in 1:M) { 
		phi[m] ~ dnorm(0, sigma[m]^(-2))
		sigma[m] ~ dunif(0,10) 
		realclus[m] <- step(sum(memb[,m])-1) 
	}
	
	# Truncated DP (alpha, G0)
	alpha ~ dgamma(5.157, 4.536)
	V[M] <- 1 
	p[1] <- V[1]
 	for (m in 1:M-1){
	    V[m] ~ dbeta(1, alpha)
		p[m+1] <- V[m+1]*(1-V[m])*p[m]/V[m]
	}
	
	beta0 ~ dnorm(0,0.01)
	for (b in 1:16) {beta[b] ~ dnorm(0,0.01)}
	gamma ~ dnorm(0.5, 0.01)
	for (d in 1:21) {delta[d] ~ dnorm(0,0.01)}
	scale ~ dnorm(1,0.01)

	cut[1] <- 0
	cut.delta ~ dexp(1)
	cut[2] <- cut[1] + cut.delta
	
	ARerr <- sum(P1[])/sum(P2[])
	
	TMatch <- sum(Match[,])/(N*7)
		
	# total clusters at each MCMC iteration			
	K <- sum(realclus[])
}

